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Continuing our investigation into the numerical properties of the Hierarchical Reference Theory, 
we study the square well fluid of range A from slightly above unity up to 3.6. After briefly touching 
upon the core condition and the related decoupling assumption necessary for numerical calculations, 
we shed some light on the way an inappropriate choice of the boundary condition imposed at high 
^■s^ ■ density may adversely affect the numerical results; we also discuss the problem of the partial differ- 

ential equation becoming stiff for close-to-critical and sub-critical temperatures. While agreement 
■ of the theory's predictions with simulational and purely theoretical studies of the square well system 

is generally satisfactory for A > 2, the combination of stiffness and the closure chosen is found to 
£jT), render the critical point numerically inaccessible in the current formulation of the theory for most 

of the systems with narrower wells. The mechanism responsible for some deficiencies is illuminated 
' at least partially and allows us to conclude that the specific difficulties encountered for square wells 

are not likely to resurface for continuous potentials. 
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I. INTRODUCTION 

r- 

In a large part of the density-temperature plane, integral equation theories are.-a reliable tool for studying thermo- 
dynamic and structural properties of, among others, simple one-component fluidsc!; unfortunately, in the vicinity of a 
liquid -vapor critical point, integral equations are haunted by a host of difficulties, leading to a variety of shortcomings 
such as incorrect and non-matching branches of the binodal, classical values at best for the critical exponents, or 
other deviations from the correct behavior at the critical singularity!! Asymptotically close to the critical point, on 
the other hand, renormalization group (rg) theory is the instrument of choice for describing the fluid; in general, 
however, RG approaches do not allow one to derive non-universal quantities from microscopic information only, i. e. 
. from knowledge of the forces acting between the fluid's particles alone. One of the theories devised to bridge the 
conceptual gap betweeji-ihese complementary approaches is the Hierarchical Reference Theory (hrt) first put forward 
by Parola and ReattoaEij: In this theory the introduction of a cut-off wavenumber Q inspired by momentum space 
O ■ rg theory and, for every value of Q, of a renormalized potential u™(r) means that only non-critical systems have to 
\ be considered at any stage of the calculation; consequently, integral equations may successfully be applied to every 
system with Q > 0, and critical behavior characterized by non-classical critical exponents is recovered only in the 
limit Q -> 0. n 

While applicability jpfJiRT to a number of interesting systems ,r*aaging from a lattice gas or Ising-modeO to various 
one-component fluidsETEl even including three-body interactional3'0, internal degrees of freedomEj, or non-hard-core 
reference systemscj, was demonstrated early on, the main focus of research on hrt has since shifted to the richer 
phase behavior of binary systemsli3lL2l. Nevertheless, in the light of hrt's high promise and low penetration into the 
liquid physics community, further study and critical assessment of this theory seem worthwhile, even and foremost in 
the case of simple one-component fluids: indeed, it is in this comparatively simple setting that we may gain important 
insights into the numerical side of the theory, and barring special mechanisms relevant to some specific model system 
only, any problems uncovered here must be expected to haunt more advanced applications of hrt, too. InjXmr-H'ork 
we have found it convenient to restrict ourselves even further, implementing hrt in its usual formulatioratJ'Ej for 
purely pairwise additive interactions via a potential v(r) obtained from the superposition of an infinitely repulsive hard 
sphere serving as reference system, u rof (r) = w hs Xr), and a predominantly attractive tail w(r), w(0) < 0. Here we have 
made use of the notation introduced previouslyO: superscripts always denote the system a quantitiy refers to (here, 
"ref" and "hs" for the reference system and hard spheres, respectively; similarly, "(Q)" for the system with cut-off 

(cf. sub-section III B of ref. |l|) 



O 



Q), and a tilde indicates Fourier transformation. In the present contribution we apply our recent re-implcmcntation 
of the theoryla to one of the simplest potentials exhibiting phase separation, viz. the square well potential -y sw [- e Ao-] 
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Considering density-independent potentials only and chosing the hard core diameter a and the well's depth e as units 
of length and energy, respectively, the attractive well's range A is the sole remaining parameter; in this report we will 
study values of A from slightly above unity up to 3.6. 

With just one parameter, viz. A, to vary, square wells obviously make for a convenient test case of hrt and, indeed, 
of liquid state theories in general; consequently, a great wasy simulational and theoretical efforts have been directed 
at this system, and studies of its phase behavior aboundElTLJ. But square wells are also of interest in their own right, 
serving as — albeit somewhat crude — models-of— a wide variety of physical systems including, e. g., 
H 2 , CO2, CH 4 , C2H6, n-pentane and n-butaneE^TEJ while current interest in this potential derives-jroa. 
finding that square wells capture the essential features of the interactions found in colloidal systemstJO. Yet anothp. 
motivation for this first application of hrt to square wells comes from a recent, very accurate simulation studytil 
confirming and quantifying the presence iru-the. system with A = 1.5 of the Yang- Yang (YY) anomaly expected and 
experimentally found for asymmetric fluidscdH. 

Due to the extensive amount of data available in the literature the more recent of which will shortly be presented 
later on, and in view of some of the limitations of hrt in its current formulation we cannot expect to gain new insight 
into the system at hand with a level of precision comparable to that of the more sophisticated simulation schemes. 
Instead, in the present contribution our focus of interest lies on some aspects of hrt's numerical side, specifically on 
those that are sensitive to the potential's range: indeed, as stated already in ref. |l^, for a potential as pronouncedly 
short-ranged as square wells some of the numerical problems should show up much more prominently than in other 
systems like, e. g., the hard-core Yukawa fluid previously consideredlij where they are, of course, in principle still 
present but do not manifest themselves as clearly. 

In accordance with the preceding remarks, another reason why application of hrt to square wells might be worth- 
while lies in the closure underlying seemingly all applications so far of hrt to simple one-component fluids with hard 
sphere reference part: As the usual formulation of hrt in these cases relies on an ansatz for the two-particle direct 
correlation function 02(7) very much in the spirit of Stall's Lowest-Order ^-ordered Approximation£3£3 (loga) or 
the equivalent Optimized Random-Phase Approximation^ (orpa) by Andersen and Chandler, the direct correlation 
function can never extend to larger r values than the potential v(r) itself. In particular, for the square well potential 
v sw[-e,X,a]^ we necessarily have cg(r) = for r > \a so that all moments of C2(fL i. e. J K3 d 3 rc2(r) r n , n > 0, 
exist, which is obviously at variance with the correct behavior near the critical pedatpj: furthermore, at intermediate 
Q the direct correlatkupJunction can hardly be considered satisfactory, especialljii3'E2l close to r = A a. While some 
earlier publications&E^TCJ already blamed unsatisfactory aspects of hrt results on this inadequacy of the closure, 
square wells should bring out related problems of hrt with the usual LOGA/ORPA-style closure much more clearly, 
and the numerical procedure's response should provide us with a signature to be looked out for in other systems, 
too; also, even within the LOGA/ORPA-style approximation the implementation of the core condition via approximate 
ordinary differential equations (odes) for the relevant expansion coefficients easily shown to be inadequate for very 
short-ranged potentialsEj casts some doubt on the range of A values amenable to an hrt treatment in the current 
formulation of the theory Determination of the admissible A-range, on the other__haud, is particularly interesting in 
the light of refs. ^7] and ^ as well as in view of the global renormaliaation scheme&LJ originally developed by White 
and co-workers as an extension of Wilson's phase-space cell methodca to the liquid state; it is only by combining tests 
internal to the theory and comparison with data available by other means that we are able to answer this question. 

In this contribution, after a sketchy presentation of the underlying theory itself (section [Hj) and the implementation 
used (section III), in section [V we turn to the results of applying hrt to square well systems of variable range. After 
a short summary of the critical point's location as obtain ed from simulation-based and other purely theoretical studies 
of square wells for various values of A (sub-section |IV A ) we first look into the core condition's implementation, which 
provides us with a first hint regarding the range of A valu es acc essible to hrt in its current formulation and once more 
highlights the decoupling assumption's role (sub-section IV B). The latter is also implicated in the correct choice of 
the boundary condi tion imposed at high density as discussed, alongside the boundary condition's location's effect, 
in sub-section IV C. — A particularly grave aspect of hrt's numerical side i s the stiffness of the partial differential 



equation (pde) for close-to-critical and sub-critical temperatures (sub-section IV D ), the v estiges of which are evident 
in the results obtained for quasi-continually varying A as presented in sub-section IV E. A short summary of our 
findings and conclusions ends our contribution (section ^) . 
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II. THE THEORY 



The definite resource on hrt ia the rcviewfl by the theory's original authors, summarizing its formalism as developed 
in a series of earlier publicationaJT3 so that. we here present only an overview of the theory used and of its formulation, 
recapitulating, some of our earlier findingslij; the notation employed here of course co-incides with that of our preceding 
contributiontj. 

As mentioned already in section |, hrt's mainstay is the implementation of the suppression of long-wavelength 
fluctuations characteristic of RG methods by means of a renormalized potential . Thus, rather than directly going 
from a reference fluid the properties of which are assumed known to the fully interacting, system, i. e. from pair 
potential v ie{ (r) to v(r) = v rci (r) + w(r), HRT proceeds via a succession of rather artificialtj intermediate potentials 
v^'(r): For every value of the cut-off wave number Q, v^' is given by 



«W)(r) = v rcf (r) +«?W)(r) 
w(k) : k > Q 
: k < Q 



w^(k) 



^sw[- e, A.er] g sin A a k—\ a k cos A a k 

where the last line specializes to the square well potential of eq. (Q); obviously, v Tef and v are recovered in the limits 
Q — > oo and Q — * 0, 

w (°°)(r) = v Icl (r) = v hs (r) 
u(°)(r) = v(r) = v sw (r) , 

allowing HRX,to gradually turn on fluctuations of ever increasing wavelength by lowering Q from oo to zero 
(numericallyta, from Qoo to Qo); as mentioned before, criticality (together with non-classical critical exponents) 
and phase separation (with isotherms rigorously flat in the two phase region) are obtained in the limit Q — > 0. In this 
procedure it is essential to maintain the differential picture implied by RG theory and to make sure that the transition 
from Q to infinitcsimally smaller cut-off Q — dQ is continuous even in the limit Q — > 0. The latter requirement ne- 
cessitates replacing the usual free energy and two-particle direct correlation function c^\r) of the hypothetical 
system with cut-off Q and potential v^(r), the "Q-system", by suitably modified quantities, viz. 

V ~ V 2 



^(O)-0W)(O)J +f(0(O) -^(0)) 

C«)(r) = 4 Q) (r) + 4>{r) - 0»)(r) 
(f> = -p w (3=l/k B T, 

where g is the number density of the system at hand; the higher order correlation functions c$\ri, . . . ,r n ), n > 3, 
are free from such problems. (Note that all the direct correlation functions including C^(r) are taken to include the 
ideal gas termsa.) 

With this set of quantities continuous even in the limit Q — ► 0, viz. and the c$\ n > 3, hrt is derived 

as a non-terminating hierarchy of coupled odes at fixed density g, calculating the properties of the Q-system by 
treating the system at infinitesimally higher cut-off Q + dQ as a reference system; of these equations, usually only the 
evolution equation for viz. 

dQ v V ) 47T» Jn l 1 C(Q)(Q)J ' { ) 



as well as the important compressibility sum-rule 



(3) 



valid for any cut-off Q directly enter practical calculations. 

When combined with a closure on the two-particle level, eqs. (||) and (||) define a pde in the (Q, £>)-plane; it is this 
PDE that we will concern ourselves with in the remainder of this text. Said closure, reminiscent of LOGA/ORPA but 
adding one free parameter to allow imposing thermodynamic consistency as embodied in eq. (||), is given, just as in 
our earlier contributionE3, by 
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C«)(r, g) = 0(r, <?) + ^(g) «o(r, e) + l^ Q \r, g) 

JC(Q)(r,g) = G^(r,g) + ^(r,g) (4) 
G (Q) {r,g) =Zn=iln Q \Q)u n (r,g), 

where we have generalized to density-dependent potentials. Basis function Uo(r, g) is chosen to coincide with@ 
w(r, g)/w(0, g), and the higher basis functions u n (r, g), n > 1, vanish outside the core; for our specific choice of basis 
functions see appendix B of ref. [l9|. In order to ensure that both the core condition, i. e. g(r, g) — for r < a(g) 
where g(r, g) is the pair distribution function, and sum-rule (^ are met it is necessary to choose the correct set of 
expansion coefficients 7n (f?), n > 0, at every cut-off Q and for every density g; assuming their validity for Q — oo 
and adopting the short hand notations 



(here, tp is an arbitrary function of k and g), both relations can be combined to 

E liQ) [uj{k,g) {u n (k,g) -u o (k,g)u n (0,g)) ,g\ 



= a W{g)±W{u j {k,g)u Q {k,g),g] (5) 

, Q 2 MQ,Q) uj(Q,g) ■ > , . 

c«3)(Q, e )(c(«)(Q, e )-0(Q,e)) ' J ~ ' 

this set of equations must, of course, be truncated to a finite number 1 + N CC of basis functions, and it is also necessary 
to neglect non-local contributions to dX^'[ip(k, g), g]/dQ to allow convenient evaluation at arbitrary Q. Both of these 
approximations have been discussed at length in our previous contributions, and while the value of N cc was found to 
strongly influence the quality of the results obtained, determination of the "fiP\g) from eq. (^) and said approximation 
for the slowly converging Z-integrals' Q-dependence leads to systematic deficiencies at small r in g(r) as determined 
from the Ornstein-Zcrnikc relation. — Unfortunately, for numerical reasons£3 it is necessary to also adopt the so-called 
decoupling assumptioncl, viz. aW)(g) = 0; as can easily be seen, this is not only mathematically incompatible with 
thermodynamic consistency but even suffices to decouple the pde implied by eqs. (||) and (||) to a set of unrelated 
ODEs at fixed density only lacking thermodynamic consistency and thus unable to predict clear phase bouiidarieslij; 
furthermore, we cannot rule out that decoupling may have a significant influence on the solution generatedlij, which 
is particularly troublesome as the much longer range of wo(r) oc (f>(r) when compared with the other basis functions 
was originally invoked as justification for setting a^\g) = 0: for square wells, this assumption is certainly even less 
justified than for the rather long-ranged hard-core Yukawa system (z = 1.8/ct) considered in ref. 

Returning to the PDE, for the numerical implementation's benefit we, too, adopted a re-formulation in terms of an 
auxiliary function f(Q, g) simply related to the modified free energy's derivative with respect to Q. The details of 
the procedure leading to a pde of the form 



^f(Q,Q) = d 00 [f,Q,Q] + d 01 [f,Q,g}JLf(Q,g) + d 02 [f,Q,g]^f(Q,g), 



9 ttn ^ = d 00 [f,Q,g] + d 01 [f,Q,ci 9 *fn ~\^j~ifn .1 a 2 
f(Q,g)u 2 Q (Q,g) = ln(l d(Q)(Qe)J , tiQ){Q . B y 

and the coefficient functions d ai themselves can be found in appendix A of ref. |l9|, q. v. ref. 

The above formulation (||) of the problem, of course still coupled to the odes implementing the core condition, 

obviously has to be amended by initial and boundary conditions; while the former easily follow from ^l?^ = 0, 
n > 0, (which is sufficient to also determine f{Q oc, (?)), choice of appropriate boundary conditions is slightly more 
complicated: if, as is the case in most of the calculations reported here (exceptions see below), the low density 
boundary is located at g m i n = 0, we can make use of the divergence of the ideal gas term —1/g in to derive not 
only f(Q, 0) = but also df(Q, 0)/dg = which alone is, in principle, sufficient to uniquely determine the solution up 
to arbitrarily high density; for computational reasons, however, it is preferable to instead only impose vanishing / at 
gmin and to supply an approximate condition for calculating / at p ma x- Among the candidates for the constraint to be 
imposed upon the solution at £ ma x in addition to the core condition there are two we should mention here: Starting 
with ref. ^, the so-called ORPA-condition, viz. 7o'^(f?max) = 0, has been used extensively; it should, however, ha-noted 
that this condition is incompatible with both thermodynamic consistency and the decoupling assumptionEj. An 
alternative first considered in our previous reporto is the decoupling assumption (g niax ) = itself; of course, this 
condition is still incompatible with the compressibility sum-rule (|^) but this is less of a problem at a boundary where 
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partial derivatives with respect to g cannot be evaluated anyway. Another option (not pursued in this contribution) 
is to give up the core condition altogether, retaining only the lowest basis function uq in the closure and thus 
effectively replacing eq. (|5|) by eq. (|3|); this has the added advantage of mathematical consistency while still retaining 
the structure of a PDE so important for obtaining clear phase boundariestj, v. s. 

It is one_of hrt's main achievements to allow calculating a fluid's binodal (coinciding with the spinodal in three 
dimensionsE3) without resorting to Maxwell const ructionsEj, for subcritical temperatures yielding isotherms rigorously 
flat in density intervals the boundaries of which are readily identified with the coexisting densities g v and gi. Thus, 
as the isothermal compressibility kt of the fully interacting system, readily found to be proportional to exp(/ — 
(0(O)//C(°))) - 1 (cf. appendix A of ref. diver ges in the two-phase region, so must the auxiliary function f(Q, g) 
in the limit Q — > 0. As a direct consequence of this, the transition from the modified free energy (g) to f(Q, g) is 
not only computationally convenient but also allows, us to follow the isothermal compressibility's build-up much more 
easily; even more importantly, a simple analysislljtj of the behavior of the PDE's coefficients for large f(Q, g) readily 
characterizes the PDE as stiff: for any density g £ [g v , gi] and close to Q = 0, the true solution f(Q, g) oscillates rapidly 
on a Q-scale of the order of exp(— /), with both an upper bound on the oscillations' amplitudes and /'s average slope 
growing like 1/Q — needless to say that this behavior cannot be reproduced numerically (v. i. sub-section IV D| ; q. v. 
ref. |l9|). Note, however, that it is not an artifact of the c&. writing of the pde in the form (JsJ) but rather a problem 
inherent to hrt itself in a formulation based upon eq. (]3)E3. 



III. NUMERICAL PROCEDURE 

The numerical study of HRT for square well systems of varying range parameter A in section FV| has only become 
feasible due to our recent re-implementatiorie£l of this theory, discussed at length in refs. [l9] and |20| ; we will make use 
of results obtained with this program exclusively. From a practical point of view, our software provides a means of 
solving a finite-difference approximation to the PDE (|^) in an iterated full-approximation scheme, imposing boundary 
conditions at densities g min and f5 max as well as initial conditions at Q — Qoo, generating a solution for Q as low 
as Qo while ensuring numerical soundness of every step by employing a number of criteria. The pivotal parameter 
governing all of the numerics is a small quantity denoted e# characteristic of the maximum admissible relative error 
introduced in a single step in the —Q direction; due to the paramount importance of derivatives with respect to g, e# 
is strictly related to the coarseness of the density grid. 

The only exception to the general strategy of ensuring a numerical quality of e# at every step in the calculation is 
the choice of step sizes AQ in the —Q direction, at least for sub-critical and close-to-critical temperatures: indeed, 
in that part of the (Q, p)-plane where the divergence of the isothermal compressibility builds up, the pde's stiffness 
(v. s.) renders fixed-precision arithmetic and relative errors bounded by e# incompatible. Consequently, for the 
calculations r y£g>orted below we resort to step sizes AQ pre-determined in a way analogous to that employed in earlier 
applicationsBO; still, monitoring and assessing suitable components of the solution vector in terms of e# as described 
in sub-section III E of ref. [l9] may yield a wealth of information on the numerical process and its evolution. 

Most of the calculations reported here have been performed on an equispaced density grid of N e = 100 density 
intervals spanning the range from g m - ln = to £> max = 1/c 3 , corresponding to a value of e# = 10~ 2 ; N cc was usually set 
to 7; and the pre-determined step sizes started from AQ — — IQ~ 2 /a at Qoo — 80/ct, plunging to a mere —5 • 10~ 6 /<7 
when approaching Qq = 10~ 4 /er. — When locating the binodal via the divergence of the isothermal compressibility 
Kt we did not require an actual overflow to occur but instead looked for a Kx-ratio at neighbouring densities exceeding 
10 4 , which is a rather reliable indicator for the binodal's location as kt typically jumps by far less than two or by 
at least some twenty orders of magnitude within one Ag; the reported values for g v and gi are the mid-points of the 
density intervals so found. In principle this allows us to locate the coexisting densities and the critical temperature 
and density to arbitrary precision, even though the computational cost rises sharply with falling e#. 



IV. APPLICATION TO SQUARE WELLS 

As mentioned before, much of the motivation for applying hrt in the formulation outlined in section [n] to the simple 
square well model potential is based upon various observations indicating possible limitations of this apnroximate 
formulation of HRT for short-ranged potentials. A case in point is the recent work of Caccamo et al£3 entirely 
devoted to several thermodynamically consistent theories' ability to deal with narrow hard-core Yukawa systems; sure 
enough, in the case of hrt the shortcomings of the LOGA/ORPA-style closure (||) and, presumably, of the accompanying 
decoupling assumption underlying the core condition's implementation via eq. (|^) were manifest already in refs. || and 
|l3| and recently confirmed by usEj, q. v. ref. 
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Of course, any of the problems discussed below only relate to hrt when implemented along the lines of sections [0] and 
[II and not to hrt proper; however, for reasons discussed in ref. [l^, alternative formulations almost certainly render 
the numerics far more demanding and open up a whole new suite of problems regarding the numerical implementation's 
soundness, especially when performing Fourier transformations of cut-off affected functionsEl 

In the following sub-sections we will complement the discussion of ref. ^| by further investigation into the numerical 
nature of hrt; before that, however, it seems pertinent to re-iterate some of the points raised in that publication as 
far as they concern the reasoning to be put forward in the following. In particular, according to section IV of ref. [l^, 
for the numerical results to be meaningful the coexisting densities g v and gi must maintain a separation of at least 
several density grid spacings Ag from the boundaries at g m - ln and g max ; consequently, j3 should never exceed some 
maximum value, f3 < /3 max , and for the systems considered here and in ref. [19| and for the typical choices for g m i n 
and g max the binodal's proximity to the low density boundary renders /3 max largely density grid- and e#-independent. 
— Not to be confused with /3 max is the lowest temperature ks / /3 ma x,# numerically accessible to the program with 
pre-determined step sizes: this is the temperature below which the program of section III never reaches Q ks Q or 
produces abnormal results; note that /3 ma x,# may be larger or smaller than /3 max , depending on the chosen combination 
of physical potential, approximations in the formulation used (the boundary conditions in particular), and the choice 
of parameters affecting the numerical work. 

Regarding the implementation of the core condition as sketched in section |n 
that a minimum of N cc = 7 basis functions in addition to uq were necessary for acceptable results despite residual 
defects of g(r) close to the origi n; a sh ort discussion of the core condition's slightly different role for square wells will 
be presented below (sub-section [V B[ ). 

The critical density g c predicted by hrt, it should be noted, is virtually always in reasonable agreement with 
literature data as shortly presented in sub-section IV A; indeed, HRT is even able to reproduce the_marked rise in 
g c predicted by refs. 



the main conclusion of ref. 19 



21, E7j, and E9| for A — > 1+ as opposed to the rigorously constant value in ref. 



Due to the 

satisfactory g c -values obtained numerically we will henceforth exclude g c from the discussion; for a demonstration of 
both gcS insensitivi ty to variation ofparameters of the numerical procedure and the quantitative agreement with the 
data of sub-section IV A see figures [j] and |[ 

In this context it may be of interest that the hrt estimate for the critical density presents no difficulties for 
the hard-core Yukawa fluid considered in ref. [H^ cither, nor is there any mention of such difficulties in any of the 
other publications on this topic that we are aware of; indeed, the theory's numerical problems primarily lie in the 
solution's small-Q behavior for close-to-critical and sub-critical temperatures on the one hand and the use of mutually 
incompatible assumptions prompted by the need to employ decoupling without giving up thermodynamic consistency 
on the other hand. Both of these aspects of hrt pertain to different parts of the (Q, p)-plane, located close to the high 
density boundary for the role mathematical inconsistencies play and at not t oo l arge Q and g ~ g c for the pathological 
behavior related to coexistence; they will be discussed in sub-sections IV C and IV D| , respectively, and their vestiges 
will also be seen in the results of appl ying h rt in the formulation of section [H to square wells of quasi-continually 
varying range A e]l,3.6] in sub-section [VE. 



A. Non-HRT results for square wells 

For comparison purposes we have collected in tables | and [n| the critical temperatures of various square well systems 
as obtained from simulations (table |) or by purely theoretical means (table |0j) ; the data included have been published 
within the last decade. 

Of the simulation based results included in table [jj only those of ref. ^6] for A G {1.25,1.375,1.5,1.75,2} have 
been obtained by molecular dynamics (MD); most of the other simulation studies rely on one or the other variant 
of the Monte Carlo (MC) method: Among these, the Gibbs ensemble MC (GEMC) calculations of ref. [2l] set out to 
determine critical exponents, (3 in particular; that work's finding of /3 <~ 1/2 for A = 2 as opposed to the expected 
(3 ~ 1/3 found for A up to 1.75 prompted re-examination of the square well fluid with A = 2 by GEMC augmented 
by finite-size scaling (FSS) techniques^, refuting the mean field value for the effective exponent. — Especially in the 
critical regime, grand canonical MC (GCMC) simulations incorporating histogram reweighting and FSS offer some 
advantage over GEMC due to the latter's restriction to fixed temperature; such an approach has been applied to 
square wells with A = 1.5 and 3 in ref. a more elaborate GCMC scheme-Hot biased towards the Ising universality 
class and taking into account the YY anomaly has recently been applied toEil A = 1.5, v. s. Yet another method goes 
under the name of thermodynamic- or temperature-and-density-scaling MC (TDSMC); it was applied to the case of 
A = 1.5 and analyzed in terms of an effective Hamiltonian in refs. ^4] and |2^. — Also included in table Q are the results 
of ref. |7[ employing an MC scheme modified to take advantage of a speed-up possible by combining simulation data 
with an analytical ansatz for the chemical potential; the efficiency of this approach originally devised to study phase 
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separation allows a large number of systems to be considered. (The error bounds given for these "modified MC" 
results in table [| have been obtained from the different results displayed in ref. ^ for different parameter settings.) 

The theoretical predictions for the critical temperature listed in table |l| comprise a second-order analytic pertur- 
bation theorytJ (APT2) applicable tqJ < A < 2 and claimed accurate for A > 1.4 as well as the hard-sphere van 
der Waals (HSvdW) equation of stateES. In addition, though not listed in table [□], we have utilized the non-square- 
well-specific Okumura-Yonezawa (OY) estimates for (3 C , primarily as a starting value when looking for the critical 
temperature in our hrt calculations; for w sw [-^^] ( the OY prediction is k B T c /e = 0.203 (2tt/3) A 3 - 0.273. 



B. The core condition 



Ever since application of hrt to continuous fluids started, the implementation of the core condition has been a majpt. 
issue, probably motivating adoption of the closure (Q) and variants thereof for non-hard-sphere rfifeE ence . systemsEil 
despite its known deficiencies in the first place; indeed, it is no coincidence that several studies&OEJtlJ primarily 
concerned with the RG aspect of the theory chose to completely eliminate the core condition. When applying hrt as 
a regular liquid state theory, on the othe^Jiand, this is not an option: too great is the effect this may have on both 
correlation functions and phase behaviorEl From table III where we compile the critical temperature T c = X/ks Pc 
for various square well potentials as functions of the number N cc + 1 of basis functions in the closure (^), just as in 
ref. [w] we find virtually constant critical temperatures for 1 < iV cc < 4; on the other hand, the amount of variation 
seen upon further increasing N cc strongly depends on A, which immediately carries over to the pair distribution 
function g^°'(r, g) and its compatibility with the core condition: For A = 3, the longest ranged potential considered 



in table III, g^°\r, g) = 0, r < a, holds reasonably well except very close to r = even for N cc — 1; when increasing 
the number of basis functions all the way to N cc = 10, the pair distribution function has to be corrected for very small 
r only, yielding a \g^°\r, g)\ that remains bounded by some 10 -2 of the contact value g^°\cr+, g) for all r < <r; the 
corresponding small change in g^°\r, g) and C^°\r, g) is reflected in the near-constant predictions for /3 C evident 



from table III. Similarly, for A = 1.5 and A = 2 and within the -/V cc -range considered, the implementation of the 
core condition does not convincingly improve except for supercritical temperatures and intermediate densities; this 
time, however, the pair distribution functions remain far from compatible with the core condition even for A^ cc = 10, 
and neither /3 C nor g^°'(r,g) itself nor, for that matter, the final values of the loga/orpa expansion coefficients 
7n (ff) indicate that the expansion (Q) for C^\k, g) might be close to convergence. But if the quality of g^°\r, g) 
improves only little if at all, the remaining deficiencies are probably to be blamed on the approximation for the poorly 
convergent integrals' derivative with respect to Q mentioned earlier (cf. eq. (12) of ref. |l9|) rather than on an insufficient 
number of basis functions; on the other hand, even though the decoupling assumption cannot directly affect the pair 
distribution function's compliance with the core condition, the approximation of neglecting the non-local term in 
[tp(k,g),g] jdQ is on the same level as that of setting a^\g) = 0, as was stressed by the authors of ref. || 
upon jointly introducing these two assumptions. Thus, combining the above findings regarding the core condition 
with the analogous analysis of sub-section IV of ref. [jl] and with that contribution's investigation into the decoupling 
assumption's possible effect (cf. fig. 2 of ref. |l9|) we are led to the conclusion that decoupling poses certainly no less 
a problem here than for the hard-core Yukawa potential studied there. 



C. High-density boundary condition 



Numerically, there are two ways for the implementation of section III to fail to reach Q = Qq, both, of course, easily 



detected by the "monitoring" variant of our code ( cf. sub -section III E of ref. |1S|): due to the solution's pathological 



behavior wherever f(Q,g) is large (cf. sub-section 1VD), or because of inappropriate boundary conditions at high 
density. As for the latter — an issue intimately linked to the decoupling assumption — , the immediate reason for the 
program's failure is a near-discontinuity in the numerical solution close to the boundary: For the moment setting aside 
the decoupling assumption and other approximations, in the application of hrt with the closure (|J) at any point (Q, g) 

in the interior of the pde's domain the core condition uniquely determines the j^\g), n > 1, for given 7o^(f?); this 
expansion coefficient is then determined by imposing thermodynamic consistency as embodied in the compressibility 
sum-rule (^). At a boundary, i. e. for g G {g m m> ftnax}i however, the second density derivative cannot be evaluated 
reliably so that some other condition must be imposed; in the calculations reported here (with the obvious exception 
of those for fig. |4|) we always choose g m i n = so that the divergence of the ideal gas term in provides f(Q, q) = 
as a convenient and unproblematic boundary condition. For g = g maK , on the other hand, we are in principle free 
to use any suitable approximation for the structural and thermodynamic properties of the Q-system and to calculate 
f(Q, £max) from said approximation, thereby providing the necessary boundary condition for the PDE (|); but for 
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practical reasons it is desirable to use the same LOGA/ORPA-form for the Q-system's direct correlation function at g max 

as in the rest of the problem's domain so that, in particular, the loga/orpa prescription 7g^' (g max ) = is a natural 
choice of boundary condition. In general, however, due to the PDE's diffusion-like character any condition imposed 
at £> max that is incompatible with the solution for g < g max by necessity induces a corresponding near-discontinuity 
in f(Q, g) close to the boundary; within the framework of a finite difference scheme this is reflected in a mismatch 
°f /(Q,f?max) and the solution at densities close by, i. e. f(Q,Q max — «Ap) for small i > 1, and the mismatch's 
severity may serve as a direct measure for the inappropriateness of the boundary condition at gmax in relation to the 
approximations applied at densities in ]p m i n , gmax[- 

On the other hand, the numerics become intractable unless we adopt the decoupling assumption, and the only way to 
consistently r»se (g) = without abandoning the core condition is to decouple the HRT-PDE to a set of odes at fixed 
density only£3; this, unfortunately, removes-all traces of thermodynamic consistency from the equations and thereby 
precludes obtaining clear phase boundaries^. It is therefore necessary to restrict decoupling to the implementation of 
the core condition only while retaining the structure of a pde together with the compressibility sum rule (|^) despite 
the latter's incompatibility with decoupling. Thus, for p min < g < g max , both C^(0,g) = —d 2 ((3A^/V) /dg 2 
and a^\g) = are used for different parts of the problem; at p ma x, however, again any approximation allowing 
calculation of f(Q, Q max ) may be used, so that it is tempting to once again resort to the LOGA/ORPA-condition of 

vanishing 7g^' ) (g m ax) or variants thereof. 

But due to the decoupling assumption's possibly large effect, any boundary condition that does not incorporate 

a^Hfi'max) = — and bear in mind that 7o < ^(Pmax) and cx^(g max ) cannot both vanish at the same time for generic 
cut-off Q — will once again incur a fatally large mismatch; if, however, we must resort to decoupling anyway, it seems 
preferable to consistently apply it for the boundary condition rather than to inconsistently combine it with a condition 
alien to the theory; also, though the mismatches' magnitudes from imposing a^\g max ) — alone or from mixing it 
with the loga/orpa condition 7^ (f> max ) = generally do not differ much as long as the pde's stiffness does not 
play a role (e. g., for A = 1.5, as long as we restrict ourselves to Q ~ 8/a or higher, or to f3 <C f3 c ), the former approach 
fares better than the other one more often than not. It is only in this sense, i. e. presupposing a LOG a/ ORPA-like 
ansatz even at g max and application of decoupling in the implementation of the core condition according to eq. (@) at 
all g, that the results are largely independent of the choice of boundary condition as claimed, for (3 < (3 C , in ref.^. 

In the numerical work we find that such a mismatch is present whenever the calculation proceeds via mathematically 
inconsistent or conflicting approximations; in the case of square wells with their comparatively short potential range, 
however, the problems are much more severe than in othe r syst ems so that /? maXi # is rather small and even drops below 
P c for most of the A interval from 1 to 2 (cf. sub-section IVE ). Restricting ourselves to (5 < /? max ,# and Q = Qo, the 
mismatch is typically reflected in an increase by one order of magnitude in the three-point finite-difference estimate 
of, e. g., \d 2 f(Qo, g)/dg 2 \ right at the boundary over the near-constant values at slightly lower densities; apart from 
a positive correlation with e#, the mismatch's severity is qualitatively unaffected by a change in parameters of the 
numerical procedure or the choice and location of the boundary condition (with the above provisions). 

Another effect worth mentioning in connection with the boundaries is the influence their locations, viz. g m i n and 
£>max, may have. The basic mechanism and its implications for the coexisting densities were already mentioned in the 
opening remarks of this section; here we only want to point out that the non-criticality enforced by the boundary 
conditions not only may unduely distort the binodal predicted by hrt as demonstrated in fig. [l], very small £> m ax 
may also allow one to reach Q = Qo at higher /3, thus effectively raising /3 ma x,# while lowering /3 m ax- — Sometimes, 
however, the expectation of the binodal keeping a separation from the boundary of several Ag at least does not hold, 
and a preposterous two-phase region appears very close to g max or, very rarely, close to g mul ; e. g., for A = 1.88 and 
f3 = 0.392/e the equations can be solved all the way down to Q = Qo — 10~ 4 /ct, predicting an unrealistic two-phase 
region extending from 0.845 (5)/<r 3 to 0.995(5) /a 3 . This behavior turns out to come from the interplay of the mismatch 
at £ max and the numerical treatment of the stiffness of the pde. 



D. The region of large f(Q, g) 

For subcritical temperatures the HRT-PDE's true solution's erratic behavior in that part of the (Q, g)-plane where 
f(Q, g) is large and the isothermal compressibility's divergence is built up (cf. section O, q. v. ref. |20| ) obviously eludes 
reliable numerical realization; in particular, while e# still characterizes the level of accuracy in auxiliary calculations, 
the same can no longer be true for the accuracy of the pde's discretization as this would require step sizes AQ so 
small as to cause floating point underflow upon evaluating, e. g., Q — (Q — AQ), thus rendering finite differences 
numerically insignificant. 

Consequently, in this respect we have to give up our strategy of controlling the numerical procedure so as to locally 
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ensure a quality-.of e# at least, turning to pre-determined step sizeala AQ in addition to fixed Ag, to which similar 
concerns applyEj; on such a coarse mesh of (Q, g)-points underlying the finite difference scheme, however, the true 
solution cannot even be represented adequately, and the numerical approximation for f(Q,g) obtained from the pde's 
discretization with these far too large step sizes cannot be trusted to faithfully represent even the average behavior 
of/(Q,g). 

This inadequacy of the step sizes is reflected in various peculiarities of the solution vector obtained in the numerical 
procedure; indeed, when monitoring the evolution of f(Q,g) and the core condition coefficients Jn(g), our code 
readily detects the plummeting step sizes necessary and signals the incompatibility of the behavior seen with the 
assumption of smoothness underlying finite difference schemes. Another telltale sign is iterated corrector steps' 
failure to converge when f(Q, g) is large: even though implicit schemes like the one we employe!] are the standard 
treatment for stiff systems, the rapid growth of the oscillations' amplitudes renders the finite difference equations 
themselves unstable under iteration; only when resigning on any control of the numerical error and refraining from 
iterations of the corrector step do the step sizes AQ chosen allow one to force advancing Q all the way to Qo in 
remarkably many cases. Also, comparison of f{Qo,g) as obtained with different sets of step sizes AQ reveals that, 
for g v < q < gi, the evolution of f(Q,g) seen numerically is, driven by the number and size of Q-steps only and 
certainly does not correspond to an average over oscillations^!; the same mechanism is also responsible for a small 
AQ-dependence of the critical temperature /3 C . — By the same token, due to the efoi- and do2-terms in eq. (^|), the 
pde's stiffness and the related problems have a direct bearing on the solution outside the coexistence region even if 
the numerical predictions there turn out rather insensitive to variation of parameters of the numerical procedure; in 
particular, we expect a gradual but non-negligible distortion (in addition to the effects of numerical differentiation 
close to the near discontinuity) of the binodal, increasing with falling temperature. 



E. hrt results for square wells 



In the light of the preceding exposition as well as of the discussion in ref. [19] it may at first seem surprising that 
HRT in the formulation of section [n] has a record of being well applicable to a variety of systems (cf. section ||); also, 
as we shall see in a moment, even for square wells, a system expected to be particularly vulnerable to the problems 
just outlined, we find reasonable es timat es o f the critical points' locations for a wide range of A- values. Still, the 
mechanisms sketched in subsections IV C and [V D as well as the difficulties presented in ref. |l9| remain and manifest 
themselves numerically in a number ways. 

For a first orientation, let us look at the results summarized in figs. |^ and || where the critical temperature T c and 
density g c are shown as functions of A; the underlying calculations have been obtained with e# = 10~ 2 , imposing 
decoupling in a consistent way at g max = l/cr 3 and with A cc + 1 = 7+1 basis functions in the expansion (^) of 
the LOGA/ORPA- function Q^K With the exception of some spurious results at A ~ 1.1, whereever /3 C < /3 ma x,# the 
critical temperature in general compares quite favorably with the data of tables || and [ilj from the calculations we have 
performed for a large number of systems in the range 1 < A < 3.6 and ignoring some isolated results, a critical point-is 
found for 1.06 < A < 1.24, for 1.45 < A < 1.53, and for A > 1.939; calculations with 7V CC = 5 yield analogous resultsEi, 
with p c < /3 max ,# in a somewhat larger part of the parameter range, viz. for 1.09 < A < 1.58 and for A > 1.896, but 
will not be considered in the following in view of the considerations of sub-section IV B and of other defects that turn 
out to be larger than for N cc — 7. 

For the moment setting aside the data for A < 1.939, hrt's predictions for the critical temperature are generally 
found to be in satisfactory agreement with the /3 c (A)-curve expected from the simulation-based and theoretical results 
presented in sub-section [V A. Embedded into this regular overall behavior of (3 C as a function of A, however, we find 
a number of depressions and elevations of (3 C , some of which cannot be seen on the scale of the plot |^ but from the 
numeric results onlyc9; others, however, are so strong as to render the critical temperature a non-monotonic function 
of A, which is certainly not expected on the grounds of the literature presented in sub-section IV A, the data of 
refs. |27H29] in particular. 



In the light of sub-section IV D it is of course tempting to simply attribute this behavior to the difficulties previously 
discussed, especially since the critical point is located in the region of large f(Qo,g) by definition; the peculiar 
distribution of A-values affected, however, suggests that these problems of the numerical procedure are triggered by 
a special mechanism. Indeed, a closer look at the core condition function C^'(Q, g) for fixed density g reveals, for 
every single one of the A values implicated that we checked, that the combination of terms pertaining to w(r) or v hs (r) 
alone (of ranges A a and a, respectively) regularly and quite frequently reduces the amplitude of this function's swings 
about the ideal gas value of — 1/g; the same happens only occasionally for A-values removed from these irregularities 
so that it is, in fact, possible to quite reliably determine whether or not a given A is affected from a plot of C^(Q, g) 
for g ~ g c alone as illustrated in fig. 0. It will come as no surprise that most of the irregularities occur when A, 
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the ratio of the two characteristic lengths present in the model, is close to a simple fraction: among the shifts in T c 
most obvious are those where A is close to 2 (cf. fig. |3|), 2j, 2-|, 2^ and 2-^, and in retrospect it seems justified to 
also include the small parameter range around A = l\ in this list, v. i.\ the effect is less obvious from fig. ^ but still 
discernible at 2|, 2^ and 2|, whereas for 2^ and 2 1 it is so small as to make the plot of /3 C (A) appear smooth while 
the irregularities are still evident from the numerical values; also note that, once again, g c is hardly affected. 

All these observations seem to indicate that indeed it is the interplay of the two different length scales and the 
resulting partial oppression of a significant portion of the oscillations of C^'(Q, g) that cause the discrepancy of hrt 
and literature results for the critical temperature around certain A values even though a smoot h inte rpolation of hrt's 



predictions from A values nearby is well compatible with the data presented in sub-section [IVAJ. Even though we 
currently cannot pinpoint the precise mechanism by which this unphysical behavior of hrt arises and, in particular, 
cannot distinguish between the closure's inadequacy and the pde's stiffness as the main culprit — though the latter 
is certainly implicated to some degree — , two conclusions may be drawn quite safely: for one, as long as we stay clear 
of values of A > 2 that are close to simple fractions, or restrict ourselves to A > 2.7 where the effects are rather small, 



we can probably trust the numerical results — with the caveats of ref. |19| and sub-sections IV C and IV D — to the 
same degree of confidence as those obtained for the hard-core Yukawa system in ref. |l9|. And secondly, it is only in the 
presence of discontinuities in the potential that certain lengths feature prominently in the relevant functions' Fourier 
transforms and can so give rise to problems of the kind outlined above; consequently, as long as we confine ourselves 
to continuous w(r, g), which still includes most of the potentials popular in liquid state physics, the unphysical shifts 
in j3 c seen for certain parameter combinations are likely not an issue, whereas the same problems are expected to 
resurface, e. g., for the multi-step potential also defined in sub-section III B of ref. [l9]. 

Another lesson to be drawn from the findings presented here as well as in ref. ^| is that, as a general rule, conclusions 
should never be drawn from isolated results alone; it is only through the combination and meticulous scrutiny of a 
set of related calculations that meaningful information can be extracted from HRT calculations: due to the problems 
related to the implementation of the core condition, to the nature and location of the boundary conditions, and to 
the pde's stiffness, any single calculation must be considered as of uncertain standing. As an examples, the analogue 
of fig. |l| for A = 1.5 shows a considerably larger variation in the binodal and the critical point's location, which is 
consistent with the above conclusions regarding the reason for /3 max ,# rising above (3 C in a narrow region around this 
A value, whereas any one of the phase boundaries found in itself is a perfectly plausible candidate for the "true" HRT 
binodal. 

V. CONCLUSION 

In conjunction with the findings of ref. |l^, the discussion of section [jV] provides quite coherent a picture of hrt's 
numerical side as well as of some peculiarities encountered for square wells. Most prominently, we see a marked 
dependence of the quality of-the results on the potential's range, confirming the trend^Qf— decreasing accuracy for 
narrower potentials reportedc3 for the hard-core Yukawa fluid; it has long been acceptedtlE^O that the simplistic but 
computationally convenienO closure (|J) has a part in this, and an improved closure has recently been proposedEil. 
Still, as far as numerical application of hrt is concerned the closure cannot be discussed without reference to the 
decoupling assumption and to the approximate implementation of the core condition via odes coupled to the hrt-PDE; 
and while the former has been found problematic both for square wells (present contribution) and for the hard-core 
Yukawa fluid considered in ref. |l^ and should probably not be trusted easily for any system, the severity of the 
difficulties brought about by the simplified treatment of the core condition sensitively depends on the potential type 
and parameters chosen: for the continuous and rather long ranged Yukawa potential with z = 1.8/a, g^°'(r, g) can 
be made sufficiently smal l with in the core, and the s quar e well fluid with A = 3 fares equally well at least; from the 



discussion of sub-section IV B and the data of table [II]|, however, it becomes apparent that smaller A — we have 
looked at A = 1.5 and A = 2 in particular — incurs substantial problems, with residual defects in the pair .distribution 
functions attributed to the ill-justified approximation of neglecting some slowly converging integralsElu in the core 
condition. 



But table III demonstrates not only the A-dependence of the results' sensitivity to the number 7V CC + 1 of basis 
functions retained in the truncated eq. (||) when varying N cc in the range < A cc < 10: while the virtually constant 
critical temperature predicted for A = 3 seems trustworthy and is, indeed, well compatible with simulation results 
(cf. table |), the amount of variation in (3 C for A = 1.5 and, to a much lesser degree, for A = 2 precludes accurate 
determination of the critical temperature; this is a first indication that the theory might be able to handle square 
wells with A = 3 quite reliably whereas problems cannot be denied for A = 2, and A = 1.5 seems largely out of reach 
for hrt in the present formulation. This is confirmed by the results obtained by quasi-continuous variation of A in 
the range 1 < A < 3.6 as shown, for A cc = 7, in fig. 0: the critical point is accessible only in part of this parameter 
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range, and not only the critical temperature at fixed A but also the bound aries of the A-intervals where HRT is able 



to reach temperatures as low as T c strongly depend on iV cc (cf. sub-section [VE ) . 

Applicability of hrt to only a restricted A-range is, of course, again related to the pronounced short-rangedness 
of the square well potential; to explain it, however, we have to invoke not only the LOGA/ORPA-style closure and 
the approximate implementation of the core condition but also the other difficulties encountered in the numerical 
procedure as highlighted in this and our prec eding contribution, viz. the decoupling assumption (ref. |l9|) , inappro- 
priate boundary conditions (sub-se ction tVC ) and, most importantly, the PDE's stiffness for thermodynamic states 
of high compressibility (sub-section [V D| ). All of these are, in principle, always present to some degree in numerical 



applications of hrt; it may prove valuable that sub-sections IV C and IV D provide distinct signatures readily de- 



tected by the implementation of section III that might be looked out for in more advanced applications of the theory, 



too. — Related to these difficulties is a peculiar effect specific to square wells (sub-section EVE): close to certain 
A values, simple fractions in particular, we see shifts in the critical temperature that render hrt's predictions much 
less compatible with simulations and other theoretical descriptions of the square well fluid than would be expected 
from the results obtained for A values close by; the mechanism for triggering these local distortions of the /3 c (A)-curve 
is illuminated at least to the point of linking it to the presence of a discontinuity in the potential's perturbational 
part. All in all, the numerical evidence as well as comparison with literature data suggest that the formulation of hrt 
sketched in section || is well able to deal with square wells and to locate their critical points to reasonable accuracy 
for A <; 2 as long as certain values are avoided, or else for A > 2.7. 
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1.219(8) GEMCEiL 
1.222 TDSMCE 
1.226 TDSMCE 
1.246(5) TDSMCE 4 
1.27 MDE3 r 
1.302(8) mod. MC 


Y)f 

.2" 

,21 
,2E 

7 




1.65 


1.645(5) mod^ 


MCt 




1.75 


1.79 MBtfU n 
1.811(13) GEMCBr 




1.8 


2.062(8) mod^ 


MCfU 


2 


2.61 MDlfJ 
2.648(14) GEMC+FSS 
2.666(85) GEMC+FSS 
2.678(27) GEMC+FSS 
2.6821(8) GEMC+FSS 
2.684(51) GEMC+FSS 
2.721(89) GEMC+FSS 
2.730(14) GEMC+FSS 


22 
22 
22 
22 
22 
25 
22 
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2 764(231 
2.778(7) 


GEMCB 
mod. MC 


27 




2.2 


3.80(1) 


mod. MC 


2T 




2.4 


5.08(2) 


mod. Mi^ 


2', 




3 


9.87(1) 


GCMCO 



TABLE I. The critical temperature T c of square well systems for various values of A as predicted by simulations and 
simulation-based theoretical analyses, and the corresponding references. The acronyms used for labeling the method employed 
in obtaining these results are given in sub-section (VA of the text. 



A 


k B T c (X)/e 


method 




1.125 


0.587 


APT2 


2(1 




1.25 


0.751 
0.850 


HSvd" 
APT2 


M 

29 


i 


1.375 


0.978 
1.08 


HSvd" 
APT2 


M 

29 


l 


1.5 


1.249 
1.33 


HSvd" 
APT2 


M 

K 


i 


1.625 


1.61 


APT2 


ill 




1.75 


1.859 
1.93 


HSvd" 
APT2 


29 


i 


1.85 


2.23 


APT2 


2!) 




2 


2.506 
2.79 


HSvd" 
APT2 


2' 





TABLE II. The critical temperature T c of square well systems for various values of A as predicted by purely theoretical 
means, and the corr espon ding references. The acronyms used for labeling the method employed in obtaining these results are 
given in sub-section [VA of the text. 





k B T c (X= 1.5)/e 


k B T c {\ = 2)/e 


kg T C (A = 3)/e 




1.209437(035) 


2.660946(132) 


9.891032(298) 


1 


1.190663(034) 


2.682489(105) 


9.899937(478) 


2 


1.203326(035) 


2.686289(105) 


9.900894(478) 


3 


1.200152(035) 


2.686078(105) 


9.900894(478) 


4 


1.197136(034) 


2.685655(105) 


9.900894(478) 


5 


1.287443(040) 


2.527365(093) 


9.737080(462) 


6 


1.098329(029) 


2.742404(110) 


9.822071(471) 


7 


0.984757(047) 


2.914763(124) 


9.867502(475) 


8 


1.070878(027) 


2.744830(110) 


9.773324(466) 


9 


1.216333(036) 


2.749695(110) 


9.887510(477) 


10 


1.207583(035) 


2.937591(126) 


9.748203(464) 



TABLE III. The dependence of the critical temperature of square well systems on the number A cc + 1 of basis functions 
retained in eqs. (fl) and (o) for various values of A. For A^ cc > 0, the decoupling assumption was imposed as high density 



boundary condition, whereas the LOGA/ORPA-condition 7o ((?max) = served the same purpose for A^c = 0; other parameters 
were chosen as indicated in section fill. 
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"(0.01, 0.5) 
*(0.01, 1.0) 

111 02 CK3 0^4 
ga 3 

FIG. 1. The b inodal of the square well system with A = 3 as obtained for different values of (e#, £Wx), cf. the discussion in 
sub-section [VC. Note that for this rather long-ranged system the critical point's location is virtually unaffected by variation 
of these parameters. Also, imposing the boundary condition at Q max = 0.5/a 3 clearly induces a shift in g v to higher and, to a 
lesser degree, in qi to lower values even well above the temperature where gi gets close to £> max , which is readily interpreted as 
an effect brought about by stiffness; the results for g max = 0.5/a 3 and e# = 0.005 (not shown in the plot) do not differ much 
from those with the same £? max and e# = 0.01 except in the binodal's vapor branch's shift being somewhat smaller. 
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FIG. 2. The critical temperature T c (dots in upper panel) and critical density g c (bars in lower panel) of square well systems 
for A ranging from close to unity up to 3.6 as obtained from calculations with parameters chosen as indicated in se ction [II; 
also included are the non-HRT predictions listed in tables Q and |fl| labeled by the acronyms introduced in sub-section IV A and 
already used in those tables. The ticks on the top border of the figure's frame indicat e the A values considered; of the 200-odd 
systems we looked at, /3 max ,# exceeded /3 C only in the A ranges indicated in sub-section 

those ranges. The three boxes in the upper panel indicate the parameter ranges displayed at larger scale in fig 
panel, the bars show the coexisting densities found according to the prescriptions of section III for the highest-temperature 
sub-critical isotherm calculated in locating the critical temperature, which explains the apparent differences in g c 's accuracy; 
the smallest g c intervals shown coincide with the spacing Ag = 10~ 2 /<r 3 of the density grid. 
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FIG. 3. The critical temperature data of fig. |2| for values of the square well range parameter A close to 1.1, 1.5, and 2 at 
larger scale; the symbols coincide with those used in fig. |^. 




FIG. 4. The core-condition function C^(Q, g) for g = 0.3/u 3 , j3 = 0.2/e and for two different ranges A of the square well 
potential, on arbitrary scales; the horizontal lines correspond to the ideal gas value —1/g. Note that, for A = 3 (upper curve), 
the peak of every single one of the function's swings is partially reduced, which is the case less than half the time — and at 
rather high Q only — for A = 2.9 (lower curve). We have excluded the data for Q < 10/a so that the effects of the pde's 
stiffness are still negligible; the underlying calculations have been performed by solving the odes corresponding to consistent 
application of the decoupling assumption at the density indicated. 
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